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The quasi-harmonic (QH) approximation uses harmonic vibrational frequencies ujq^h{V) com¬ 
puted at volumes V near Vb where the Born-Oppenheimer (BO) energy Eei{V) is minimum. When 
this is used in the harmonic free energy, QH approximation gives a good zeroth order theory of 
thermal expansion and first-order theory of bulk modulus, where n*’'-order means smaller than the 
leading term by e", where e = /jUvib/J^ei or ksT/Egi, and E^i is an electronic energy scale, typically 
2 to 10 eV. Experiment often shows evidence for next order corrections. When such corrections are 
needed, anharmonic interactions must be included. The most accessible measure of anharmonicity is 
the quasiparticle (QP) energy ujq{V,T) seen experimentally by vibrational spectroscopy. However, 
this cannot just be inserted into the harmonic free energy Eh- In this paper, a free energy is found 
which corrects the double-counting of anharmonic interactions that is made when E is approximated 
by FH{ijJQ{V,T)). The term “QP thermodynamics” is used for this way of treating anharmonicity. 

It enables (n-l-l)-order corrections if QH theory is accurate to order n. This procedure is used to give 
corrections to specific heat and volume thermal expansion. The QH formulas for isothermal (Bt) 
and adiabatic (Bs) bulk moduli are clarified, and the route to higher-order corrections is indicated. 


I. INTRODUCTION 

In non-magnetic, insulating materials, thermodynamic 
behavior is controlled by vibrational excitations, which 
are often close to harmonic. There is a unique and correct 
version of harmonic theory, based on Taylor-expanding 
the Born-Oppenheimer (BO) energy £'ei({i?f}) around 
the atomic coordinates {Ri^o} of a crystal with volume 
V. The BO energy is a “ground state” property, and 
the main target of density functional theory (DFT). 
This gives normal mode eigenvectors and frequencies 
^Q,h{V), where the index Q = {Q,j} labels the states; 
Q runs over the N wave-vectors of the Brillouin zone, 
and j runs over the 3n branches. These states can be 
called “non-interacting quasiparticles.” However, in this 
paper, the term quasiparticle (QP) will be reserved for 
the vibrational resonances seen experimentally. These 
differ from harmonic eigenfrequencies because higher or¬ 
der (“anharmonic”) terms in the Taylor expansion are 
not negligible. 

In this paper, the term “harmonic” will refer to the 
unique correct harmonic limit, further specialized to the 
case when the volume is chosen to be Vo, where the BO 
energy is minimum. It is useful also to know the har¬ 
monic normal modes (and their frequencies ujq ,Hiy)) 
at other volumes; this is the “quasi-harmonic” (QH) 
theory. The harmonic approximation is a good start¬ 
ing point, successfully implemented by “a6 initio” DFT 
calculationJ3, and useful, often to good approxima¬ 
tion, for things like specific heat, Cp{T). Vibrational 
spectroscopy^ of reasonably pure crystals most often sees 
reasonably sharp Lorentzian resonances. They can be as¬ 
signed a wave-vector Q, and are expected to show a one- 
to-one correspondence with the harmonic normal modes. 
They are the QP’s of this paper. The central frequency 
LOQ (the QP frequency, or energy) is T-dependent. There 


is good evidence from theory-experiment comparison^ 
that the QH energy wq{V), evaluated at the correct 
thermally-expanded volume V(T) at higher T, does not 
reproduce well the measured thermal shifts Alcq of QP 
energies LOgiV, T). There is an anharmonic thermal shift, 
additional to and different from, the pure volume-related 
shift of QH theory, and it is often significant at higher T. 

Terminology is not universally agreed upon. Cowlej^^, 
in his seminal paper, derived the modern Matsubara per¬ 
turbation theory for anharmonic effects. He occasionally 
uses the word “quasi-harmonic” to denote what is here 
called “quasiparticle”. In recent literature, QH most of¬ 
ten denotes use of ojq^h{V) from T = 0 DFT. Occa¬ 
sionally papers about anharmonic theory still use “QH” 
and “QP” interchangeably when referring to approximate 
normal modes, ujq{V,T), which are here called QP. 

The QP relaxation rate 1/tq is the full width at half 
maximum of the spectroscopic Lorentzian line. In pure 
crystals it lies outside harmonic theory, and is also T- 
dependent. This paper is about extracting additional 
thermodynamic information from the temperature de¬ 
pendence of u}Q, ignoring 1/tq. This suffices for most 
low-order thermal corrections. I will call this “quasipar¬ 
ticle thermodynamics”. It differs from “quasi-harmonic 
thermodynamics”. Deviations from harmonic vibrations 
are responsible for thermal shifts of many properties. 
The ones of prime concern in this paper are Cp and Cy 
(constant pressure and constant volume), bulk modulus 
Bp and Bs (isothermal and adiabatic), and volume ex¬ 
pansion V{T) and a = {ljV)dVjdT (constant pressure.) 
Good general discussions are in older literature.®^ 

There are two main ideas in QP theory: (1) low- 
lying excitations correspond 1-to-l with a non-interacting 
single-particle picture; they have QP energies fujjQ{V,T) 
and mode occupancies (ng); and (2) low energy dynam¬ 
ics can be described as the dynamics, in space and time, 
of the mode occupancy. QP theory can fail in at least 
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two ways, (a) The resonance can be very non-Lorentzian 
so that UJQ is poorly defined, (b) Correlated occupancy 
{fiQfiQi) — {nQ){nQ>) may become important. In this pa¬ 
per, UQ denotes the equilibrium (Bose-Einstein) mean oc¬ 
cupancy [exp{fujjQ[y,T)/ksT) — , n-Q.H its harmonic 

or QH version, and (fig) denotes the actual occupancy 
in a general ensemble, not necessarily equilibrium. En¬ 
tropy plays a special rol^, since it is just 1/N times the 
logarithm of the number of ways of distributing {n)Nhuj 
of excitation energy into N oscillators of frequency w, 

S = kB ^[((fig) + 1) ln((ng) -hi) - (ng) In(ng)]. (1) 

Q 

The equilibrium occupancy ng is the one which max¬ 
imizes Eq.Q at fixed energy, and the thermodynamic 
entropy S(T) is given by Eq.Q with (ng) —)■ ng. When 
harmonic frequencies ujq^h are used in ng, the result 
is the harmonic entropy Sh- When the T-dependent QP 
energies are used in nq, and inserted in S(T), an accurate 
improvement of the thermodynamics is achieved. This 
will be denoted S’qp. The same replacement does not 
work for the harmonic free energy. If ujq(T) is inserted 
into Fh, the resulting F does not obey —dF/dT = Sqp. 
A “correct” QP free energy that does agree with S'qp is 
found as follows. The thermodynamic energy U{T) = 
E -|- TS" in harmonic theory is 

Uh = 'y^^ ^Q,H{nQ,H + 1 / 2 ). ( 2 ) 

Q 


When normal mode frequencies acquire an anharmonic 
correctionj _ i^q,h —>■ ^q{T), the energy acquires a 

correction,E2l 

^QP = y^^a;g(r)(ng -|- 1/2) 

Q 

-(l/2)^%g(r)-a;g.ff](ng + l/2) (3) 

Q 

This corrects for double-counting of the interaction, but 
only in leading anharmonic approximation wher^ 


a;g(T) - ujQ^H = = 



( 4 ) 

where the superscript indicates the lowest order cor¬ 
rection (second-order perturbation theory) which comes 
from third and fourth-order anharmonicity. Further de¬ 
tails are in the Appendix. Here dojq/dnqi is a T- 
independent anharmonic interaction function. In higher- 
order perturbation theory there are presumably addi¬ 
tional shifts of the type 



( 5 ) 

involving anharmonic interactions up to sixth-order in 


displacement uq. We assume these can be omitted, and 
this approximation enables the correction in Eq.([^ to be 
sufficient. Then an accurate and consistent thermody¬ 
namic free energy is Fqp = Uqp — TSqp. Notice that 
the anharmonic shift, Eq.Q, does not vanish at T = 0, 
but has a zero-point component where ng/ -1- 1/2 —>■ 1/2. 
The quasiparticle frequencies are shifted from the har¬ 
monic frequencies even at T = 0. 

QH thermodynamics is a limiting case. It uses only 
volume-dependence of a;g ij(P). The correction factor 
in Uqp, Eq. (§, vanishes, and the quasi-harmonic free 
energy is just the harmonic free energy with a volume- 
dependent harmonic frequency. It improves pure har¬ 
monic theory and gives correct lowest-order thermal 
shifts for properties such as the bulk modulus which are 
volume derivatives of F. The reason it works to lowest 
order is because dujq/dV differs from dujq^n/dV only in 
next order. QH theory is computationally accessible, but 
QP theory much less so. QP theory suffers from the dif¬ 
ficulty that the anharmonic shift is not usually measured 
except at a few temperatures. It can be numerically com¬ 
puted using DFT for the anharmonic forces. It is not yet 
computed routinely, but this is changin^A^Ii^. 

Exact theory associates vibrational resonances with 
poles of a phonon Green’s function, a correlation func¬ 
tion describing dynamics on the BO energy surface. Ex¬ 
act thermodynamics should be computed from a corre¬ 
sponding theory for the free energy. This can be com¬ 
puted perturbativeljl41ill. At high T, classical molecular 
dynamics (MD) describes dynamics non-perturbatively, 
if the BO forces are known. This is called ab initio MD, or 
AIMD. Thermodynamics gen erally then requires a tricky 
“thermodynamic integration”E2E®^. To do a correct 
non-perturbative computation at lower T requires quan¬ 
tum corrections, as in path-integral MD.^^ 

Zero-point and related thermal nuclear motions cause 
shifts and iso tope-depe ndences in many measured phys¬ 
ical propertieP"^I^2H28l Explicit formulas are given here 
for the first-order shifts of Cp, Bp, Bs, and a. If no spec¬ 
ification (adiabatic versus isothermal) is made, isother¬ 
mal is implied. The adiabatic shift can be found by ther¬ 
modynamic rules, as shown in Sec. lYg Many of these re¬ 
sults can be found in some form in the literature. There 
is a lot of correct,I^M^ plus much partially correct, as well 
as incorrect or confusing literature on this subject. There 
are semi-empirical formulas that have a long history of 
enabling useful fitting, even th ough the formulas do not 
seem to be justifiable in detail.^SEni xhe aim of this pa¬ 
per is a simplified, possibly less confusing, derivation of 
correct formulas. 

The paper is organized as follows. In Sec. |H| exam¬ 
ples of thermal shifts from experiments are given. In Sec. 
|TTTj extra complexities of non-cubic crystals, and crystals 
with internal coordinates, are discussed. In Sec. |IV[ the 
QH approximation and the QP approximation are dis¬ 
cussed. Specific heat formulas are presented, showing 
that QP theory provides a thermal correction inacces¬ 
sible in QH approximation. In Sec. two orders of 









3 


thermal correction to the volume are discussed. This 
gives Griineisen theory of thermal expansion ao plus a 
first order correction. In Sec. |VI[ the leading correction 
to the bulk modulus is derived (from QH theory). The 
Appendix reviews the microscopic theory of Eqs.(l]|4). 


II. EXPERIMENTAL EXAMPLES 


Figures illustrate the thermal shifts under dis¬ 
cussion. Fig. [ij shows that the bulk modulus has sur¬ 
prisingly large vibrational correction^^IHini, T hese have 
serious significance for geoscience, for example.l^^llll The 
leading-order bulk modulus, Bq, comes from electronic 
stiffness. The product BoVa (14 = Vo/Nn is the volume 
per atom in leading order theory) has order of magnitude 
10 eV, a characteristic electron energy. Vibrational en¬ 
ergies are two orders of magnitude smaller. I will define 
e as the dimensionless ratio of phonon to electron ener¬ 
gies. This will appear explicitly in the form Hcu/BVa or 
ksT/BVa in various results. A parameter like e « 0.01 
controls the size of the vibrational corrections under dis¬ 
cussion. Fig. [IJshows « 10% shifts, indicating that there 
can be a significant prefactor multiplying e. 



FIG. 1. Experimental bulk modulus B versus temperature. 
Isothermal data for ice Ih are from ref. [311 For NaCl, the 
low T adiabatic data are from ref. 1321 and high T data are 
from ref. M The isothermal data for NaCl are from ref. IMl 
For MgO, adiabatic low T data are from ref. |35|and high T 
data from ref. 1361 Isothermal data for MgO are from ref. 1371 
The MgO data are larger by 10 than the others, and have a 
weaker thermal shift. Many metals (Al, Cu, Ni, etc., ref. I38D 
have bulk moduli similar in magnitude to MgO, with large 
thermal shifts (similar to NaCl, for AT ~ 0 d). Covalent ma¬ 
terials, like carborP^and silicorP^, are more similar to MgO, 
showing weaker T-dependence . The dashed lines are approx¬ 
imate extrapolations suggesting that the zero-temperature, 
frozen-lattice values of B are larger than the measured zero- 
temperature values by ~ 15% (NaCl) and ~ 9% (MgO). The 
extrapolation for ice Ih is not given, and would require de¬ 
tailed calculation as explained in the text. 


The small parameter e can be estimated as e « 
ks^D/B^Va, where 0 _d is the Debye temperature. Fx- 


perimental B and 14 may be used. Rough values are 
e = 0.0045 (silicoiM, 0.0087 (MgCPU), 0.0072 (NaCP), 
and 0.042 (ice I H^^ I ^^ I). However, the parameter e for 
ice Ih is poorly defined. The number e=0.042 used the 
low T 0£) « 300K. This measures only thermally ex¬ 
cited (acoustic and librational) vibrations at T < 273K. 
The “0-H stretch” vibrations at the opposite end of the 
spectrum have hui/kB larger by 11. These modes also 
contribute to the zero-point shifts in ice. If optic modes 
are used to define e, the value of this “small parameter” 
increases to 0.5. 



FIG. 2. Experimental volume 14 per atom versus tempera¬ 
ture. Data for silicon are from ref. 1431 Data for MgO are from 
m Data for ice Ih are from ref. Data for NaCl were con¬ 
structed by integrating the polynomial expressions for linear 
thermal expansion given in ref. 1451 The dashed curve extrap¬ 
olates the NaCl V {T) using the quasi-harmonic high T slope 
(valid for T > On) dVa/dT = Sksy/Bo, with experimental 
values 7 = 1.57 and Bq = 23.7GPa, as tabulated in ref. 1461 
This suggests that the zero-point expansion AV{T = 0)/Vb 
of NaCl is about 5%. Crude extrapolations for MgO and sil¬ 
icon suggest zero-point expansions of less than 2% and 1% 
respectively. The extrapolation for ice Ih is not given. De¬ 
tailed calculations for ice show that (as is true for silicon as 
well) Griineisen parameters are negative for some modes and 
positive for others; the theoretical zero-point expansion of ice 
Ih was computed to be « 1% in ref 1471 


It is perhaps worth mentioning that the representa¬ 
tion of a physical property P as an expansion in e 
(P = Po-l-Pie-l-P 2 e^-|-...) is not forced to have universally- 
defined coefficients. Especially because the expansion is 
truncated after the Pi or possibly P 2 term, it is normal 
that the last coefficient may (or may not, depending on 
the source) contain some higher effects (Pi = Pio -l-Piie, 
for example.) The only rule is that Pnc" should contain 
nothing of lower order than e". This will be mentioned 
again in Secs. [V|and 

The volume shift^SEHUl shown in Fig. i are smaller 
in relative size. Fig. shows that a, the temperature 
derivative of the volum^^, roughly follows a harmonic 
specific heat (Cq) type of T-dependence. This is the 
result of Griineisen theorj®!. But at higher T, there is 
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a very significant thermal shift of both Cp{T) and a(T) 
away from the Cii{T) form. 

Figs. and also illustrate zero-point shifts. Mean 
square thermal lattice displacements of the atom, in 
harmonic theory, ar^ 


•?) = E 




*|Q)P ( '^Q,H + 


( 6 ) 


where {i\Q) is the component of the phonon eigenvector 
IQ) on the atom, and where ng is the equilibrium 
occupation number. The zero-temperature value (u^) ^ 
h/2Mu} is the quantum zero-point motion, which depends 
on nuclear mass, whereas the high-T value ksT/Moj"^ is 
classical and depends only on the force constant 
not on the nuclear mass M. The low T (u^) causes 
zero-point shifts of atomic volume V{T = 0) and bulk 
modulus, which differ for different isotopes. Therefore, 
the DFT (“frozen-lattice”) value Vq or Bq should differ 
from the actual value F(0) or B{0). It is interesting that 
the frozen-lattice value can sometimes be deduced from 
experiment.!^ This is because the thermal factor n-|- 1/2 
at high T becomes a;(l —l/12a:^-|-- • •) where x = ksT/huj. 
An asymptotic linear-in-T fit to n-|- 1/2 at high T passes 
through 0 at T = 0. It is of course difficult to find the 
“correct” experimental asymptote, since thermal correc¬ 
tions enter to alter it. However, the bulk modulus simpli¬ 
fies the fit if both isothermal and adiabatic versions are 
available, because each should extrapolate to the same 
T = 0 value Bq, as is shown on Fig. Curves of this 
type are in the review of Leibfried and Ludwig.!^ 
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FIG. 3. Experimental volume expansion a{T) and specific 
heat Cp{T), for crystalline argon, at P =1 atmosphere, up 
to the meltin g te mperature (82.3K). The data were compiled 
by Bodryakov^, who analyzed ~ 10 different experiments. 
Shown for comparison is the Debye model with Qd = 81.2K 
fitted to Cp data. Both experimental and Debye-model spe¬ 
cific heats were scaled, by the same factor, to lie on top of the 
volume expansion curve at lower T. 

We have already seen that the anharmonic shift 
Awg = ujQ — OJQ^H of phonon frequencies has a simi¬ 
lar form, Eq. 0 . The T-independent coupling parameter 
dujqldnQi has contributions like 14 and \Vz\^/hw, see 


the Appendix, Eq.(A3). They have order-of-magnitude 
etujjQ. For example, the term of the type |14p/Iia; in¬ 
volves a third order anharmonic coupling coefficient V 3 , of 
structure u^d^Eei/du^, and of order {u/a)^Eei. The ra¬ 
tio {u/aY of lattice displacement to interatomic distance 
is of order huj/Muj'^a? « e. Putting all these factors to¬ 
gether, we see that dujq/dnQr « etuq (or Aujq/ujq « e.) 

The smallness of the anharmonic shift is only a crude 
estimate which sometimes may fail. A failure is likely 
to cause anharmonic broadening Fg of vibrations to be 
bigger than the spacing of vibrational levels wg. In such 
cases, phonon quasiparticles are poorly defined, pertur¬ 
bative treatments may fail, and thermal shifts may not 
be described well by quasiparticle theory. 

The validity of perturbative computation beyond har¬ 
monic approximation for thermodynamic properties is 
not a closed issue. Wallac^^ summarizes evidence 
for failure of anharmonic perturbation theory to re¬ 
produce apparently reliable MD. But Boltzmann equa¬ 
tion treatments of thermal conductivity are now very 
successfuPSHSIl^ and are based on the third-order term 
in the same perturbation theory. Computations based 
on DFT anharmonic forces are becoming more common, 
and generally claim decent agreement with experiment. 
A nice example is theory and experiment for aluminum 
by Tang et oZ.!^. A thermal conductivity k > 5 W/mK 
is a good hint that phonon quasiparticles are mostly well 
behaved. This crude estimate comes from k ~ Cviji 
where the specific heat is C = Skp/^a^ is the volume 
per atom, v ~ TrWmax/a is the sound velocity, a the lattice 
constant, and i the phonon mean free path. Quasiparti¬ 
cle theory requires i large compared to a, or k large com¬ 
pared to Kmin = Cva/3. If we choose fiWmax/fcs 0_D 
to be 300K, and a to be SA, then ^ IW/mK, and 
K > hKmin should be sufficient to trust quasiparticle the¬ 
ory for most of the phonons of the material. However, 
0D may be significantly bigger or smaller than 300K, 
and the criterion could be scaled to k > (0£)/3OOK)x 5 
W/mK. 


III. NON-CUBIC CRYSTALS AND INTERNAL 
COORDINATES 

Pressure, volume, and temperature are not the only 
thermodynamic variables in crystals. One can also have 
anisotropic stress Gap and anisotropic strains Cap- Pres¬ 
sure and volume change are the traces of these tensors. 
This paper looks only at pressure and volume. The gener¬ 
alization to tensor properties complicates notations and 
results, but the principles are not changed. Consider 
hexagonal structures as a simple example of non-cubic. 
The separate a{T) and c(T) lattice parameters are rel¬ 
evant thermodynamic variables. They are also not con¬ 
sidered in this paper, only V{T) = \/Za?cl2 is consid¬ 
ered. When T changes, not only does V change, but 
also c/a. This can not be ignored, but is kept hidden 
in this paper. The volume-dependent phonon frequency 
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ujQo{V) is treated as a well-defined variable. There is a 
hidden assumption that this has been computed at var¬ 
ious volumes, and for each volume, the correct c/a ra¬ 
tio has been found and used in the phonon calculation. 
Finally, consider wurtzite crystal structure (hexagonal 
symmetry and 4 atoms per cell). There is one “internal 
coordinate” u{T), which determines the c-axis offset be¬ 
tween the cation and anion sublattices. This also cannot 
be ignored. But it is hidden by the assumption that for 
a particular choice of V, the correct u{T) as well as c/a 
have been computed, and used to find uiq^h = WQ(Vb), 
and dujQ^H/dV, etc. Cubic crystals can also have inter¬ 
nal coordinates not fixed by symmetry, which need to be 
treated the same way. 


IV. QUASIPARTICLE THERMODYNAMICS 

Even when harmonic approximation is seriously per¬ 
turbed by anharmonic effects, there may still be phonon 
quasiparticles, with effective interactions not too strong. 
Thermodynamics is then approximated by using the 
QP frequencies ujq{V,T) in the non-interacting entropy 
formulcP, 

'S'qp = fcs X! + 1) “ ’^ q ] ■ C^) 

Q 

Because of the V and T-dependence of the QP energy, 
^Qp has altered V and T-derivative s which g ive correc¬ 
tions in thermodynamic calculations^ESESES] gQj._ 

responding free energy is 

Fqp = EeiCF) + Fvib.H(C, T) + AFah, (8) 


Pvib.H 


fcsT^ln 

Q 


2 sinh 


( hu;Q{V,T) \ 

V SfcsT J 


(9) 


AFah = - -b 1/2), 


( 10 ) 


where Ag^ = a;Q(U, T)—a;Q,jy was defined in Eq.(j^. The 
part Evib.H has the standard form of the harmonic free 
energy, but here the quasiparticle frequency ujq{V,T) is 
inserted. The double-counting correction AFah is nom¬ 
inally smaller by e than the part EVib.K- This version 
of Fqp is the same as Uqp — TSqp and Eq.(j^ for Uqp. 
The quasi-harmonic formulas Sqh and Fqh are the same 
except that ujq{V,T) is replaced by u}q{V), usually cal¬ 
culated by DFT. In that case, the anharmonic term, 
Eq.(10l, vanishes. In a metal or a magnetic material. 


one should include additional terms in Fqp for thermal 
excitation of electrons or magnons. Such effects are omit¬ 
ted here. The QH procedure of using just a volume- 
dependent QP energy in the harmonic free-energy for¬ 
mula, does give correct first-order U-derivatives, but fails 


to give thermal shifts which depend on T-derivatives. In 
this sense, it can be regarded as an incomplete, rather 
than an incorrect, theory, and as a partially correct sim¬ 
plification of QP theory. It correctly contains the infor¬ 
mation available from DPT calculations of the frequency- 
spectra at different volumes. Ramirez et made a 
careful study of the accuracy of the QH approximation 
by comparison with well-converged path-integral MD for 
three phases of ice. They find generally very good agree¬ 
ment between QH and PIMD 

As an example of QP thermodynamics, consider the 
specific heat, C = TdS/dT. The free energy is not 
needed; the correct QP entropy is Eq. 0 with QP fre¬ 
quencies in the equilibrium occupation functions. 


cv = r I'A 




Q 


dnq 

dT 


( 11 ) 


Here X is pressure P or volume V. This gives 


Cx.QP = ^ 
Q 





( 12 ) 


where the subscript “H” means {dnQ/dT)p = 
[fnoQ/kBT^)‘nQ(nQ + I), obtained by differentiating uq 
by the explicit T in the Bose function, but not by the im¬ 
plicit T contained in u}q(V, T). The first term of Eq.( 12) 
is a harmonic specific heat Ch, but not the purely har¬ 
monic Co, because the frequencies ujq appearing in the 
formula are the renormalized T-dependent quasiparticle 
frequencies. The difference between Ch and Co is a gen¬ 
tle T-dependent stretching of the harmonic Co(T) curve 
along the T axis. This does not affect the high T classi¬ 
cal limit, SNnkB- A serious high-T deviation from har¬ 
monic theory (see the measurement for Ar, Fig. must 
be caused by the second term of Eq. (HI- 

In QH theory, {dujQ/dT)v = 0, so Cv,qh = Ch- Also 
in QH theory, {dojQ/dT)p = {dujQ/dV)T{dV/dT)p, so 
there is a significant QH correction to Cp. QH theory 
gives the correct difference, Cp — Cy, but it misses the 
anharmonic correction which appears in the correct QP 
theory for both Cy and Cp. Computational evidence^ 
shows that QH theory shifts Cp(T) away from the har¬ 
monic value Cg, but that experiment exhibits different 
shifts.l^ 


V. VOLUME EXPANSION 


The aim is to get corrections to one higher order than 
the standard Griineisen theorjB^. The method is to use 
Eq.(|8l for the free energy, calculate P{V,T) = —dF/dV, 
and then find the volume V (T) at which the pressure is 
zero. It is convenient to have a notation for the dimen¬ 
sionless volume expansion, Q 

C = {V- Vo)/Vo 


(13) 
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where Vq is as usual the volume that minimizes the frozen 
lattice (Born-Oppenheimer) energy. For results to order 
it is necessary to know the DFT frozen lattice energy 
Eel to third order in C, 


EeliV) = EeliVo) + IboVoC^ + IboVoC^ + ■■■ ( 14 ) 

2 0 


Here Bq = Vo{cPEei/dV‘^)vo is the order e° electronic 
contribution to the bulk modulus, and Bq is the third 
derivative, VQ{d?Ee\/dV^)vo- Bq is similar in order of 
magnitude to Bq- The notation B is used because the 
notation B' means dB/dP — —1 — B/B in equation-of- 
state theory.[22llll Normally Bq < 0 is found. Crystals are 
softer when expanded and stiffer when compressed. From 
the volume derivative, we get an “equation of state,” 



1 

2 


Boe-H 

Q 


OEh dujQ 
dujQ dV 


dV 


ipQ + 1/2) + A 


( 2 ) dug dujQ 

^ du}Q dV 


(15) 


Making the substitutions dFn/dujQ = h{nQ + 1/2) and 
dnq/dojQ = —{h/kBT)nQ{nQ + 1), and setting P = 0, 
this becomes 




1 / 2 ) 


(16) 


Here the definitions have been introduced 


V dujQ 


(17) 


F aAg> 

WQ dV 


(18) 


where 7 q is the “mode Griineisen parameter,” and Sq is 
the analogous volume derivative of the correction ujq — 
UJQH- Therefore, 6q = jq - {uqh/ ujQ)"fQH, where 


V 


IQH = - 


^QH 


dujQH 

dV 


(19) 


of the volume (C = (C — Vq)/Vq) is 



= h ( 21 ) 


where Cqh{T) = {hujQ/V)dnQ/dT is the specific heat 
per harmonic mode. These are the famous Griineisen re¬ 
lations. Griineisen’s paper^^ of 1912 and 1918 were a 
remarkable advance, simultaneous with the first true un¬ 
derstanding of crystals that came from Rutherford and 
Bohr, von Lane and the Braggs, Born and von Karman, 
Eucken, and Debye. Geophysicists and others like to 
define an average Griineisen parameter 7 and to write 
Eq.(15) as P — Pei+^Uvih/V where Bvib is the harmonic 
vibrational energy, Eq.Q. This is called the “Griineisen 
equation of state.” It omits the higher-order corrections 
which are now to be discussed. 


B. Quasi-harmonic theory 


The need for corrections is evident from Fig. show¬ 
ing large high-T deviations of thermal expansion relative 
to specific heat. If Eq.(21| were correct, this would re¬ 


quire unphysically large T-dependence of Griineisen con¬ 
stants. To correctly find the volume shift to next order, 
it is necessary to solve Eq.(16) self-consistently. The job 


is complicated by the fact that 7 q as defined in Eq.(|l 
and jQH as defined in Eq.(19l, are volume dependent. 


An interesting example is the computation by Skel¬ 
ton et al^ of a{T) for PbS, PbSe, and PbTe. Above 
the Debye temperature, a shows strong T-dependence, 
similar to argon in Fig. These computations had no 
anharmonic corrections. This shows that a good quasi¬ 
harmonic theory does in fact have important corrections 
beyond the lowest order Griineisen theory. The volume- 
dependent electron and vibrational free energy in Eq.® 
are included, omitting the correction term involving Ag . 
The resulting Fqh{V,T) is minimized at fixed T, giving 
Vqh{T). This is equivalent to a self-consistent solution 


of Eq.|l 6 [), omitting the terms involving Sq and 


A 


( 2 ) 


A. Lowest-order (Griineisen) theory 


C. Full second-order theory 


( 2 ) 

Griineisen parameters 7 q are of order 1, while Ag is 
a small anharmonic correction of order ewg. Therefore 
(5 q is of order ejg. All the terms on the right of Eq.(16) 
except the first term involving jq are higher order in e. 
Therefore, the leading-order relation for the thermal shift 


It is rather messy to do the full solution to second or¬ 
der. To simplify things, consider just the high-T limit, 
where nq + 1/2 is « kBT/hujq and nq{nq + 1 ) is 
« (kBT/huiq)^. The quantum corrections are factors 
[1 ± {fujjjkBT)'^(12 + T^order +•••]. Then, to order e^. 
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Eq.(16l becomes 


c - 

2Bo^° 

ksT 


BoVoil + Co) 


EhQ(^)-T 


A( 2 ) 

^QO 

2 t^Q 0 


7Q0 (22) 


The subscript O’s indicate that quantities are all (except 
the first appearance of jq) evaluated at the frozen-lattice 
T = 0 minimum volume Vq. The lowest-order volume 
expansion, Coj at high T, from Eq.(20l or from Eq.(22) is 

keT 


Co = 


VoBo 


E 

Q 


IQO- 


The definition of 7 qo is 


7Q0 = 


a^Qo 


du>, 


QH 


dV 


(23) 


(24) 


Vo 


and the phonon frequencies loqq are similarly the T = 0 
harmonic frozen lattice values. 


We now need to expand the fully anharmonic 7 q 
around 7 qoj to first order in C- This is done in two stages. 
First we expand around ^qh, 


iQ{y) = 


V 9{ujqh + ^Q^) 
U^QH + Ag) W 



1qh{V) + 5qo 


(25) 


Subscripts 0 indicate sufficient accuracy for a first-order 
result. Next, we expand 'fQniV) around the volume Vq 
where 7 qo is defined. To do this, we need to know the 
harmonic frequencies to second order in C around Vq, 


^qh{V) = ujQo [l - 7QoC - 7 QoC^/2 H-] , (26) 


7Q0 = 


t /2 

WQO 


^ dV^ 


Vo 


(27) 


The notation 7 q used here (Eq|^ is not the same as 
^{djQ/dV) = 7Q + 7 q + 7 Q. Volume dependence of 
7 q(V) has often been neglected. If the mode Griineisen 
parameter were independent of volume, one could inte¬ 
grate to find ujq{V) = ujq{Vo){V/Vo)'^<^ . As observed 
previousljEHIIIl^ there is no justification for this except 
unwarranted optimism. After some algebra, the relation 
between 7 q// and 7 qo is 


1qh{V) = 7Q0 + [7Q0 + 7qo + 7Qo] C- (28) 


Combining the two stages, the result is 


Ag) 

IQ = 7Q0 + [7Q0 + 7qo + 7Qo] Co-7Q0 + <5 qo H-■ 

WQO 

(29) 

Finally, insert this into Eq.(22) and keep only hrst order 
corrections. The result is 


ksT V 


BoVo 


Q 


(75o + 7Qo)Co + ^- 


A( 2 ) 

0 


2 wqo 


IQO 


(30) 


The corresponding high T formula for the volume ther¬ 
mal expansion coefficient is 


a = ao + AaQH + AaAH- (31) 


The leading term, oq, is just the Griineisen formula eval¬ 
uated with frozen-lattice parameters. Its high-T form is 


ao 


ks 


(32) 


The high-T quasiharmonic correction is explicitly linear 
in T, 


AaQH = 2Tao ^( 7^0 + 7Qo) j - ^al 

(33) 

These are smaller than the leading term oq by a factor 
like Tao e. The high-T anharmonic correction is 


AaAH 


ks 




‘Vqo 



(34) 


This is also smaller than oq by one power of e. The 
anharmonic factors 6 — (A/w )7 vary linearly with T at 
high T. This is why the factors of 1/2 multiplying <5 qo 
( 2 ") 

and (AQ//a;Qo)7Qo in Eq.(30) disappear in Eq.(34| after 
taking the temperature derivative. 


The QH calculations of Karki et for MgO show 
that QH corrections, Eq.(33), can cause a large effect, 
even exceeding the experimental linear rise in a. The 
calculations of Mounet and MarzarP^ also show a sig¬ 
nificant QH linear increase of a(T) in diamond, but 
less than the shift observed experimentally by Slack and 
BartrairP^. These results indicate that the anharmonic 
part of Eq.(311 is as important as the QH part. A path 
integral Monte Carlo study by Herrero and Ramired^ 
confirms this. 





























VI. THERMAL CORRECTION TO BULK 
MODULUS 


The literature about B{T) is large because of its im¬ 
portance in geoscience. The bulk modulus is the simplest 
and most accessible part of the elastic constant matrix 
Cij , all components of which show related zero-point and 
thermal alteration. This paper focuses on B for simplic¬ 
ity, but generaliza tion to the full elasticity tensor is not 
hard! 6 | 7 | 38 | 39 |5ll6M9| The bulk modulus is dominated by 
the large electronic contribution Bq. Corrections to first 
and second order in e can be found from the quasiparticle 
free energy, Eq.® , and pressure, Eq.( 15). The first-order 


shift of the isothermal bulk modulus uses only the first 
two correction terms in Eq.(15), 


P(U, T) = -BoC - ^Boe - E + 1/2) (35) 


Q 


There is no anharmonic contribution to i? = —VdPjdV 
in first order. The simulations by Ramirez et alW^ (using 
a fluctuation formulcP for B{T) ) confirm the accuracy 
of quasiharmonic theory for phases of ice. For cases like 
ice, where volume shifts are relatively large, it is insuf¬ 
ficient to compute only low-order derivatives of energy 
and vibrational frequency. But QH theory, with 7 q(E) 
computed separately for different volumes along the QH 
V (r)-curve, has been shown to work.l^ 


Taking the volume derivative of Eq. (351, 


Bt,qh{T) = ^{Bq + BqC) 
ro 


-Pn 


Q 


hiOQ 

BoV 


( driQ 

I dT 


7q + 


nq 


7Q 


(36) 


The last term uses the identity d{u}Q^Q/V)/dV = 
ojq'^qIV'^. The next to last term uses the fact that 
VdnnJdV equals T{dnQ/dT)Yi'yQ. The first term of 
Eq.( |36[ ) is the purely electronic term, Pei- To first order 
in e it can be written as Po-l-(Po+Po)Co- The high T ver¬ 
sion of Co is given in Eq.([^, and the general expression 
is the same as the Griineisen version, Eq.([20|, except fre¬ 
quencies and derivatives are evaluated at (V^T) = (Vq, 0 ). 
Then Eq.([36|) becomes 




= E 


/ hujQQ \ 

\BoVo) 


_j. f dnqo 


nqo 



IQO 


\ dT 


Bo 

Bo 


7qo 


7Q0 


(?7) 


This equation is contained in somewhat hidden form in 
Leibfried and Ludwig. Born and Huan^ and Wallac^ 
also give this result, except altered because frequencies 
and derivatives are evaluated at {V,T). The paper of 
Karch et al.^^ gives an alternate derivation. Many in¬ 


correctly simplihed versions exist. 

The para met er e is not truly small for ice I h. For this 
reason, Eq. [^does not work particularly welpTEol pjp 
rect computation and minimization of the QH free en¬ 
ergy (Eq. I without the last term) may work. This 
has done used for many years, even in cases where 
the shifts are small enough that Eq.(37l should be 
adequate.^EHHmiZI] in cases, like ice Ih, where e is 


too large to use Eq. (371, there is no guarantee that truly 


anharmonic terms of order and higher are not as im¬ 
portant as QH terms found by direct minimization. 

It is important to distinguish between adiabatic (Ps) 
and isothermal (Pt) condition^^. The definitions are 

Bt = -V{dP/dV)T = V{d^F/dV^)T (38) 


Bs = -V{dP/dV)s = V{d^U/dV'^)s 


(39) 


where U and F are the internal energy and Helmholz 
free energy respectively. Thermodynamics gives exact 

rclationd *^‘7l23i2^38i63i73ff 7r 


I38l63l73irf5l 


Bs Cp To^BtV Ta^BsV 

jD ^ ^ '> 

Dj' L/y Lyy Op 

where Cy/V is the heat capacity per volume. The prod¬ 
uct aT is of order e, and CvT/BtV is also of order e, 
so the shift (Ps — Bt)/B is po sitive and order e. The 
full tensor version also is availabl^TElHZIl The vibrational 
corrections 5s (adiabatic) and 5t (isothermal) are both 
first order in e, and they differ from each other in the 
same order. The leading order value of To^BtV/Cv is 
sufficient for correcting isothermal to adiabatic. Using 
Eq.(21) and the harmonic specihc heat, the result is 


Bs — Bt = 


T 


J2 q ^qo 


(dnQo\ _ 
[ dT 


V 




dT 


(41) 


Figure shows approximate high T slopes (dB/dT) 
of both Bs and Pt for NaCl. In the high T limit 
where hujQ{dnQ/dT) —>• kp, Eq.(41) reduces to d{Bs — 
BT)/dT = Skp^'^/Va, where 7 = z^^q/^N. The slopes 
shown in Fig. 1 then require 7 « 1.5, in good agreement 
with other estimates of 7 for NaCl. 


Appendix A 

This appendix tries to illuminate the quasiparticle 
thermodynamics of Eqs.(l|4) by using anharmonic ther- 

Ac 


mal perturbation theory. According to Cowlej®, the vi¬ 
brational thermal Green’s function matrix, in the basis 
of harmonic eigenstates |A) with frequencies oj\,h, is 


(G 


-1 


)aA' — ~ ^^)^AA' + — fT^A'] 


(Al) 
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The eigenvalues uP' of the matrix G~^ + uj‘^1, with the 
imaginary part T omitted, are denoted They are 
the renormalized squared normal mode frequencies. If 
anharmonicity is weak, then in leading approximation 
these eigenvalues are + 2uj\A\\. At the same 

level of approximation, ujx = ujx^h + Cowley gives 

an explicit formula from lowest-order perturbation the- 
ory, for the anharmonic shift = Ag = ujq — ujq^h- 

The normal mode index A is now replaced by Q = (Q, j). 
Cowley’s formula can be written in the form 



dujQ 

dnQ> 


94 

-V^^\QQ,Q'Q') 

Q" 


1 

-h 

Wqh + Wq! -\- Wq 


1 1 

- 1 - 

WQ" -I- WQ/ - WQ WQ// - WQ/ -I- WQ 

1 

^Q" - ^Q' - 


(A3) 


This is an explicit form for Eq.Q. Here 1^(3) yG) 
are third and fourth derivatives of the BO potential 
taken around the periodic sites of the crystal of volume 
Vo, and the frequencies wq and the occupation num¬ 
ber uq! use anharmonic renormalization (computed self- 
consistently.) 

Cowley also derives the anharmonic free energy at the 
same level of perturbation theory. His answer can be 
written as 


harmonic entropy. 


Sh.o = ks + 1) “ InriQ] (A6) 

Q 

Consider what happens if the “quasiparticle entropy” 
is constructed by replacing the harmonic frequencies in 
Eq. (A6 1 by the anharmonic frequencies wq Tay¬ 

lor expanding to first order in Aq, the answer is 


AS = Sqp - Sh,o =-hY, 


9nQ .{2) 
dT ■ 


(A7) 


This is the same as the entropy dF/dT from Eq. ( [A4t . 
The factor 1/2 in A4 disappears when using Eq.(A2) 


while differentiating Eq.(A4| for AF, because duQldugi 
is symmetric in Q and Q'. An alternate derivation using 
a variational principle is given in ref. [S^ This suggests 
that the use of QP energies in the harmonic entropy for¬ 
mula may be valid somewhat beyond low-order pertur¬ 
bation theory. 

Consider then what happens if the same substitu¬ 
tion ujQ^H —>■ WQ is done in the harmonic free energy 
Fh = £q^Q.h(»^Q + 1/2) - TSh- The answer is 

AFh — AQ^(nQ-|-l/2)—TAS”. This differs from the 


correct anharmonic free energy, Eq. (A41, by not having 


the correct factor of 1/2. This is a proof of the double¬ 
counting correction that was added to the internal energy 
in Eq.([^. The correct formula Eq.(A4) does differ from 
the QP theory of Eqs.([I]|^ by a small T-independent 
term Fj^q. 
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F = Fh.o + ^Y + V2) + Fao (A4) 

Q 


QQ'Q" 

1 3 


WQ" + Wq/ -I- Wq Wq// -|- Wq/ - Wq 


(A5) 


where Fp.o is the free energy of non-interacting (har¬ 
monic) quasiparticles. Now find the corresponding en¬ 
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